Evolutionary dynamics of mutants that modify population structure

Natural selection is usually studied between mutants that differ in reproductive rate, but are subject to the same population structure. Here we explore how natural selection acts on mutants that have the same reproductive rate, but different population structures. In our framework, population structure is given by a graph that specifies where offspring can disperse. The invading mutant disperses offspring on a different graph than the resident wild-type. We find that more densely connected dispersal graphs tend to increase the invader’s fixation probability, but the exact relationship between structure and fixation probability is subtle. We present three main results. First, we prove that if both invader and resident are on complete dispersal graphs, then removing a single edge in the invader’s dispersal graph reduces its fixation probability. Second, we show that for certain island models higher invader’s connectivity increases its fixation probability, but the magnitude of the effect depends on the exact layout of the connections. Third, we show that for lattices the effect of different connectivity is comparable to that of different fitness: for large population size, the invader’s fixation probability is either constant or exponentially small, depending on whether it is more or less connected than the resident.


Introduction
Evolutionary dynamics is the study of how different traits arise and disappear in a population of reproducing individuals.Each trait might confer a fitness advantage (or disadvantage) on its bearer, thus in turn altering the probability that the trait spreads through the population (an event called fixation) or disappears (extinction).Besides the fitness advantage, another important factor in determining the fate of a trait over time (its fixation or extinction) is the spatial structure of the population [1][2][3][4][5].For instance, the population might be subdivided into 'islands': an offspring of a reproducing individual then typically stays in the same island, but occasionally it migrates to some nearby island.The fixation probability of a trait then crucially depends on the dispersal pattern, that is, the migration rates among the islands.Incorporation of population structure into a model of selection dynamics substantially improves the descriptive power of the model [1,[5][6][7][8][9][10][11].
Evolutionary graph theory is a powerful framework for studying natural selection in population structures with arbitrarily complex dispersal patterns [12][13][14][15][16][17][18][19].On an evolutionary graph (network), individuals occupy the nodes (vertices), and the edges (links) specify where the offspring can migrate.Graphs can represent spatial structures, contact networks in epidemiology, social networks, and phenotypic or genotypic structures in biological populations [12,[20][21][22][23][24].The question is then: how does a graph structure affect the fixation probability of a new mutant introduced into a background population of residents?Extensive research over the past decade has produced many remarkable population structures with various desirable properties [25][26][27][28][29].As one example, consider a mutation that increases the reproduction rate of the affected individual.Population structures that increase the fixation probability of such mutations, as compared with the baseline case of unstructured (well-mixed) populations, are known as amplifiers of selection.Many amplifiers of selection are known [30][31][32][33].
In this work, we primarily consider mutations that do not change the reproductive rate of the affected individual, but rather its motility potential.That is, we consider mutants which perceive the population structure through a graph with additional edges (or with fewer edges).In nature, an altered motility potential could arise in a variety of scenarios.We give three examples.
First, consider a species occupying a region that is split by a geographical barrier into two parts.If the mutation allows the offspring to successfully cross the barrier, the mutants will perceive the population structure as being well-mixed, whereas the residents will continue perceiving it as being split into two parts (islands).This can be due to differences in migration potential between mutants and residents.
As a second example, consider structured multicellular organisms.There, cells are arranged in symmetric lattice structures known as epithelia.An epithelial tissue may be described as a two-dimensional sheet defined by vertex points representing wall junctions, one-dimensional edges representing cell walls, and two-dimensional faces representing cells.The form of this tissue network is determined by the extracellular matrix (ECM).The ECM is a network consisting of extracellular macromolecules, collagen and enzymes that provide structural and biochemical support to surrounding cells.The composition of ECM varies between multicellular structures [34][35][36][37][38]. Thus, when discussing somatic evolution in multicellular organisms, the invading genotype might differ in what network structure it is forming [36,39].In other words, each type, in the absence of the other type, forms its own and different extracellular matrix.This leads to different alignment of cells and thus a new population structure (figure 1).
Carcinoma is yet another example of how the tissue organization of the invader and resident types can differ from each other.In this case, tumour cells normally have a highly disorganized neighbourhood structure, due to the variability in cell-cell adhesion and the lack of proper epithelial programmes among tumour cells in the tumour microenvironment [40,41].Normal epithelial cells, on the other hand, are typically organized in monolayers and form geometric lattice-like patterns [42].Additionally, aggressive tumour cells, due to their high mesenchymal marker level, can show higher motility relative to normal tissue cells [43].These two facts, namely the lack of normal epithelial structure in tumours and higher motility in aggressive cancers, can lead to substantial differences in invasion pattern and potentially affect the outcome of evolutionary processes.
In order to model differences in the motility potential within the framework of evolutionary graph theory, we represent the population structure as two graphs G A , G B overlaid on top of each other on the same set of N nodes [44].The N nodes represent different sites, each occupied by a single individual.The two graphs G A , G B represent the dispersal patterns for the mutants and residents, respectively.In other words, mutant offspring migrate along the edges of G A , whereas resident offspring migrate along the edges of G B .
We study the fixation probability ρ(G A , G B ) of a single neutral mutant which appears at a random node and perceives the population structure as G A , as it attempts to invade a population of residents which perceive the population through G B .We assume that both graphs G A and G B are connected.When the two structures G A and G B coincide, the fixation probability of a neutral mutant is equal to 1/N, regardless of the graph structure.By contrast, when the structures G A and G B differ, the fixation probability can be either higher or lower than 1/N.
There is a large body of the literature on the evolution and ecology of migration and dispersal [45][46][47][48][49], especially for population structures formed by islands (also called patches, demes or metapopulations) [50][51][52][53].The present framework provides a formal way to approach the motility potential as a genotypic quality.Moreover, it enables a direct study of both simple and arbitrarily complex population structures, of any population size.In this way, the framework is a generalization of the vast literature on migration and dispersal, similarly to how evolutionary graph theory is a generalization of the vast literature on evolution and ecology in spatially structured populations [9,12].
Among the graph-theoretical approaches, other ways to model motility and dispersal have been suggested in the literature [54][55][56].They allow for the offsprings to disperse in more complex forms and reach locations that are not directly connected to the mother location.This introduces migration potential as an independent quantity relative to the proliferation potential of the types [57][58][59][60].In those cases, the motility potential is representative of a random motion, and it is typically decoupled from the reproduction events.For example, Manem et al. [58] considered a model where beside death-birth events, individual cell movements are also included.They also considered one cell type to have different motility potential than the other.It was observed that the fixation probability of the mutants is reduced when including individual cell movements.Such random motility and motion has an anti-synergistic relationship with the proliferation potential [55,57,58].
Here we show that, in contrast to random motility, enhanced structured motility generally leads to an increase in the fixation probability of the invading mutant.We do this by considering three types of population structures.First, we consider all 112 population structures on N = 6 nodes.By numerically calculating the relevant fixation probabilities, we show that increased motility potential is positively correlated with increased fixation probability.Moreover, we mathematically prove that for any population size N the Complete graph K N is 'locally optimal' in a sense described in the next section.However, we note that the effect is subtle, and we identify specific circumstances in which making mutants more (structurally) motile actually decreases their fixation probability.Second, we consider large population structures corresponding to island models, and we show that the extent to which increased motility helps the mutant fixate can vary considerably, depending on the exact layout of the extra connections.Finally, we consider large low-dimensional lattices, and we show that the effect of altered structural motility is analogous to the effect of altered reproductive rate: in the limit of large population size, the fixation probability of a mutant is either constant or exponentially small, depending on whether it is more or less motile than the residents.This latter case of low-dimensional lattices is relevant for the study of somatic evolution in epithelial tissues and carcinoma.

Standard Moran process on a graph
Within the framework of evolutionary graph theory [12], a population structure is described as a graph (network), where nodes (vertices) represent locations (sites) and the graph connectivity defines the topology and the neighbourhood.There are N nodes and each node is occupied by a single individual.Each individual is either of type A (mutant) with fitness r A , or of type B (resident) with fitness r B .The evolutionary dynamics is governed by the standard stochastic discrete-time Moran birth-death process, adapted to the population structure: at each time point, a single individual is picked for reproduction, proportionally to its fitness.This focal individual produces offspring (a copy of itself ), and the offspring then migrates and replaces a random neighbouring individual.The probability of migration from node i to node j is given by an N × N dispersal matrix M ¼ ðm i,j Þ N i,j¼1 .Thus, for undirected, unweighted graphs (which are the focus of this work), the entries m i,j of the dispersal matrix M satisfy , if nodes i and j are adjacent, 0, otherwise: (Here degðuÞ is the degree of node u, that is, the number of nodes adjacent to u.)

Moran process on two graphs
It is commonly assumed that the dispersal matrix is independent of the two types; that is, both types of individuals perceive the population through the same population structure.Following the recent work of Melissourgos et al. [44], here we study a more general case in which the dispersal pattern depends on the type of the offspring that migrates.Thus, we consider two graphs G A , G B and the corresponding . That is, any time a type A individual reproduces at a node i, the offspring replaces an individual at node j with probability m A ij .
By contrast, the offspring of a type B individual reproducing at node i migrates to node j with probability m B ij (figure 2).We assume that both graphs G A and G B are connected.
The state of the population at any given time point is described by a list n ¼ ðn 1 , . .., n N Þ of N zeros and ones, where n i = 1 denotes that node i is currently occupied by a type A individual (mutant).The model is a Markov chain with 2 N possible states.Two of the states are absorbing, and they correspond to homogeneous population consisting purely of type A individuals ðstate n 1 ¼ ð1, . .., 1ÞÞ or type B individuals ðstate n 0 ¼ ð0, . .., 0ÞÞ.

Questions and results
In this work, we study how differences in the migration and dispersal pattern G A of mutants and G B of residents influence the fate of a single random mutant which appears at a random location.As a measure of the mutant success, we use its fixation probability under neutral drift (that is, r A = r B ).We denote this quantity by ρ(G A , G B ).It is known that whenever the two types have the same dispersal pattern (G A = G B ), the fixation probability under neutral drift is equal to 1/N, regardless of the graph structure [61].Thus, the regime of neutral drift provides a clean baseline and it decouples the effect of a difference in population structure from other effects.Specifically, we study the following questions: (i) Does increased motility increase or decrease the mutant fixation probability?(ii) Can the effect be quantified for simple, biologically plausible structures, such as island models or lowdimensional lattices?
To address the first question, in section Small graphs we numerically compute the fixation probabilities ρ(G A , G B ) for all pairs G A , G B of graphs of small size.We find that, generally speaking, increased motility potential (that is, living on a graph with more edges) tends to increase the fixation probability of the mutant.In particular, we prove that the Complete graph is 'locally optimal' in the sense that if residents live on a Complete graph and mutants live on a graph that misses a single edge, then the mutant fixation probability drops below the threshold 1/N (see theorem A.1 in appendix A).However, we also identify  special cases, in which an increase in the motility potential decreases the fixation probability rather than increasing it.This suggests that for arbitrary population structures the effects of motility on the fixation probability are complex.
Given this complexity, we proceed to study pairs of specific, biologically relevant structures.Namely, in section Dense regular graphs we consider certain population structures that correspond to island models with two equal islands.We show that two such structures with the same total number of edges exhibit a substantially different behaviour in the limit N → ∞.This implies that the effect of altered motility in dense regular graphs cannot be easily quantified in terms of a single parameter (the total number of edges).Then, motivated by tissue organization in multicellular organisms, in section Lattice graphs we consider one-and two-dimensional lattices.We show that in this setting, the difference in motility can be quantified and it has analogous effect to a difference in reproductive rate: increased motility results in mutant fixation with constant probability, whereas decreased motility causes the fixation probability to be exponentially small.

Related work
The question of computing fixation probabilities for various versions of Moran processes on graphs has been studied extensively.In principle, for any population structure the fixation probability can be computed numerically by solving a system of linear equations [62].However, since the size of the system is generally exponential in the population size, this approach is practically feasible only for very small populations, or for very specific population structures [28,63,64].For large population sizes, there exist efficient approximation algorithms either in the limit of weak selection [18,29,65,66] or when the underlying graph is undirected [15,67].The two-graph setting was introduced recently by Melissourgos et al. [44], who considered it primarily from the computational perspective, and extended the above algorithmic results to a special case of mutants with substantial reproductive advantage (r A ≫ r B ) which perceive the population as a Complete graph (G A = K N ).They also considered certain game-theoretical perspective in which residents and mutants can choose their own graph structure, and they established bounds for certain specific pairs of graphs, such as the Complete graph invading the Star graph.By contrast, in this work we consider the two-graph setting from a biological perspective, and we present structural results concerning mutants with no reproductive advantage (r A = r B ) who, similarly to the residents, perceive the population structure either as an island model or as a low-dimensional lattice.On top of that, we answer a question stated in [44] related to the bestresponse dynamics in the space of all graphs.Namely, we show that while the Complete graph is locally optimal (see theorem A.1 in appendix A), it is not always the best response (figure 9).

Small graphs
In this section, we consider population structures on N labelled nodes, for small values of N. In this regime, the fixation probability ρ(G A , G B ) can be computed exactly, by numerically solving a system of 2 N linear equations following standard methods [62].
For N = 2, there is only one connected graph and, by symmetry, the fixation probability of a single type A individual is equal to 1/2.For N = 3, there are four undirected graphs: a single graph G 0 with three edges (equivalently a Complete graph, or a Cycle), and three different graphs G 1 , G 2 , G 3 with two edges each.The corresponding fixation probabilities are given in figure 3b.Note that ρ(G A , G B ) = 1/N when G A and G B are identical, but in general ρ(G A , G B ) could be both more than 1/N or less than 1/N, even when G A and G B are isomorphic (if they are not identical), see figure 3c.
For general N, there are 2 N 2 ÀN pairs of graphs on N labelled nodes.Already for N = 6 this is more than a billion pairs, hence in what follows we focus on the case when one of the graphs G A , G B is a Complete graph, denoted K N .We use a short-hand notation ρ(G) = ρ(G, K N ), for the fixation probability of a single mutant which perceives the population structure as a graph G and invades a population of residents which perceive the population structure as a Complete graph K N .Analogously, we denote by r w ðGÞ ¼ rðK N , GÞ the fixation probability of a single mutant living on a Complete graph K N and invading a population of residents which live on G. Figure 4 shows ρ(G) and r w ðGÞ for all undirected graphs on N = 6 nodes, based on the number of edges in G.

Maximal and minimal fixation probability
Among the graphs on six nodes, fixation probability ρ(G) is maximized when G is the Complete graph K 6 .Recall that ρ(K N ) = 1/N, for any integer N. In relation to this, we prove that ρ(K N ) is 'locally maximal': that is, we show that if one edge is removed from the Complete graph K N , then the resulting graph Similarly, we prove that K N is locally minimal with respect to r w ðGÞ:

Relation to the number of edges
In general, fixation probability ρ(G) tends to be higher for graphs G with more edges.However, this is only a rule of thumb.For instance, the Lollipop graph LP 6 has a relatively low fixation probability ρ(LP 6 ), given its number of edges.Here a Lollipop graph, denoted LP N , consists of a Complete graph on N − 1 nodes and a single extra edge connecting the last node.Moreover, adding edges to a graph G to produce a graph G 0 sometimes does not increase the fixation probability but rather decreases it: this is illustrated by the Pan graph P 6 and the Treetop graph TT 6 for which we have ρ(P 6 ) > 0.071 and ρ(TT 6 ) < 0.065.Here a Pan graph, denoted P N , consists of a cycle on N − 1 nodes and a single extra edge connecting the last node.In a Treetop graph, denoted TT N , the node with degree 3 is further connected to all other nodes.(We note that for the intermediate graph I 6 obtained by adding only one edge to P 6 , the fixation probability ρ(I 6 ) ≈ 0.070 is in between those two values.)

Regular graphs
Recall that a graph is regular if all its nodes have the same degree (that is, the same number of neighbours).We make two observations.First, recall that when both the mutants and the residents perceive the population structure as the same regular graph R N , the isothermal theorem of [12] states that the mutant fixation probability is equal to (1 − 1/r)/(1 − 1/r N ), where r is the mutant fitness advantage.That is, in terms of the fixation probability, all regular graphs are indistinguishable from each other.One could hope that the same holds for mutants which live on regular graphs and invade the residents on the Complete graph.That is, one could hope that for any two regular graphs However, figure 4 shows that this generalization does not hold: for two different regular graphs R N and R 0 N (even with the same degree) the fixation probabilities ρ(R N , K N ) and ρ(R N 0 , K N ) are generally different, as witnessed by the two 3-regular graphs with N = 6 nodes and 9 edges.Second, figure 4 shows that, given a fixed number of edges, the fixation probability ρ(G) is higher for regular (or almost regular) graphs as compared with non-regular graphs.For instance, among the connected graphs with six edges the fixation probability is maximized by the Cycle graph C 6 which is regular.Similarly, among the graphs with five edges it is maximized by the Line graph L 6 which is almost regular.Here a Cycle graph, denoted C N , is the connected graph where each node is connected to two neighbours, and a Line graph, denoted L N , is the Cycle graph with one edge missing.However, we prove that the Line graph L N generally does not maximize the fixation probability among the connected graphs with N − 1 edges (socalled trees): in particular, direct computation for N = 8 shows that the graph G 8 consisting of three paths of lengths 2, 2 and 3 meeting at a single node satisfies ρ(G 8 ) > 0.0098 > 0.0095 > ρ(L 8 ).

Dense regular graphs
As suggested by figure 4, regular graphs G have high fixation probability ρ(G), compared with other graphs with the same number of edges.Here we consider certain simple regular graphs that contain approximately half of the total possible number of edges.We show that for some such graphs, the fixation probability is comparable to that of a Complete graph, whereas for other graphs it is substantially smaller.Thus, the isothermal theorem [12] does not generalize to the setting with two different graphs.
Given a population size N (with N even), let B N = K N/2,N/2 be a (complete) Bipartite graph with equal parts N/2, N/2 and let T N be a Two-clique graph obtained by adding N/2 matching edges to a union of two disjoint Complete graphs of size N/2 each, see figure 5a.Note that both B N and T N have precisely (1/4)N 2 edges, which is roughly half of the edges of K N .(The graph K N has (1/2)N(N − 1) ≈ (1/2)N 2 edges.)Also, note that both B N and T N represent populations  royalsocietypublishing.org/journal/rsif J. R. Soc.Interface 20: 20230355 subdivided into two large islands: in the case of B N , the offspring always migrates to the opposite island, whereas in the case of T N , the offspring mostly stays in the same island and it migrates only rarely (namely with probability of the order of 1/N).We prove that ρ(B N ) > 0.58/N (see theorem B.1 in appendix B).Since ρ(K N ) = 1/N, this implies that missing roughly half of the edges only reduces the fixation probability by a constant factor, independent of the population size N.In fact, numerical computation shows that N • ρ(B N ) ≈ 0.82, whereas for the Two-clique graph we observe N • ρ(T N ) → 0, see figure 5b.
The intuition for this distinction is as follows.On both graphs, the state of the system at any given time point is completely described by the frequencies N L ∈ [0, N] and N R ∈ [0, N] of mutants in the left and the right half.On B N , the two frequencies remain roughly equal throughout the process (N L ≈ N R ): indeed, once say N L ≫ N R , more mutant offspring is produced on the left and they migrate to the right, thereby helping balance the numbers again.By contrast, on T N the mutants migrate rarely, thus the lineage produced by the initial mutant remains trapped in one half for substantial amount of time.Throughout that time, the mutants are 'blocking' each other from spreading more than they would block each other if they were split evenly between the two halves: indeed, with all mutants in one half, the probability that a reproducing mutant replaces another mutant (thus not increasing the size of the mutant subpopulation) is twice as large, as compared with the situation where the mutants are evenly split.For small mutant subpopulations, this effect is non-negligible and it causes the fixation probability ρ(B N ) to decay faster than inversely proportionally to N.
Regarding r w , we observe N Á r w ðB N Þ % 1:11 and the data on N Á r w ðT N Þ is unclear, see figure 5c.While we believe that N Á r w ðT N Þ % 1:4, it is conceivable that N Á r w ðT N Þ tends to a larger constant or even grows unbounded.The intuition is that when mutants live on a Complete graph K N , the offspring is equally likely to migrate to any location.By randomness, the condition N L ≈ N R is thus maintained throughout most of the early stages of the process.Therefore, as with ρ(B N ), both r w ðB N Þ and r w ðT N Þ are inversely proportional to N. To sum up, the graphs B N and T N show a considerably different behaviour in terms of ρ but a qualitatively comparable behaviour in terms of r w .

Lattice graphs
Here we study sparse regular graphs, specifically lattice graphs.Lattices exist in any number of dimensions.We focus on one-and two-dimensional lattices, since those are biologically relevant and amenable to our computational techniques.For each dimension, we study the effect of increased or decreased connectivity (degree) of the lattice on the fixation probability of an invading mutant.

One-dimensional lattices
In one dimension, we consider circulation graphs Cir d N (already studied in this context from a different point of view, see [44]).For a fixed even integer d, a d-Circulation graph, denoted Cir d N , consists of N nodes arranged in a cycle, where each node is connected to d other nodes, namely the next d/2 nodes and the previous d/2 nodes in the cyclic order, see figure 6a.
To shorten the notation, we denote by N , Cir d2 N Þ the fixation probability of a mutant living on a one-dimensional lattice Cir d1 N with degree d 1 versus a population of residents living on a one-dimensional lattice Cir d2 N with degree d 2 .Note that when When the degrees d 1 , d 2 of the mutant and resident graph differ, the fixation probability crucially depends on which of the two degrees is larger.When the mutant graph has a lower connectivity (d 1 < d 2 ) then r 1D N ðd 1 , d 2 Þ tends to 0 exponentially quickly as N → ∞, see figure 6b.By contrast, when the mutant graph has a higher connectivity (d 1 > d 2 ) then r 1D N ðd 1 , d 2 Þ tends to a positive constant c that depends on d 1 and d 2 , see figure 6c.Specifically, for large N we observe that r 1D N ð4, 2Þ % 0:16, r 1D N ð6, 2Þ % 0:17 and r 1D N ð6, 4Þ % 0:09.
Those results are in agreement with bounds 0:11 r 1D N ð4, 2Þ 0:25 that we prove analytically by a stochastic domination argument (see theorem C.1 in appendix C).The intuition behind the argument is that once the mutants form a contiguous block of a large size, the block is more likely to expand rather than to diminish at both interfaces.Indeed, royalsocietypublishing.org/journal/rsif J. R. Soc.Interface 20: 20230355 the probability of gaining the boundary node is the same as losing the (other) boundary node but, on top of that, mutants could skip the boundary node, invade the interior of the resident territory and only after that gain the skipped node.This event has a non-negligible probability of happening, hence there is a positive bias favouring the spread of mutants.For a formal proof, see theorem C.1 in appendix C.

Two-dimensional lattices
In two dimensions, we consider graphs drawn on a square lattice with periodic boundary condition.For instance, by connecting each node to its four closest nodes (Von Neumann neighbourhood), we obtain a graph Sq The results are analogous to the case of one-dimensional lattices.When the mutants live on a less connected lattice, their fixation probability tends to 0 exponentially quickly.By contrast, when they live on a more densely connected lattice, their fixation probability tends to a constant as the population size N tends to infinity (figure 7).

Effective fitness compared with complete graphs.
The behaviour of the fixation probability for pairs of lowdimensional lattices is reminiscent of the behaviour of the fixation probability ρ(K N ; r) of a single mutant with relative reproductive rate r ≠ 1 in a well-mixed population of N − 1 other residents.In that setting, we have ρ(K N ; r) = (1 − 1/r)/ (1 − 1/r N ).For any fixed r ≠ 1, this formula exhibits one of two possible behaviours in the limit N → ∞.When r < 1 then ρ(K N ; r) decays approximately as 1/r N .By contrast, when r > 1 then it tends to a positive constant 1 − 1/r.(When r = 1 we have ρ(K N ; r) = 1/N by symmetry.)This is a similar dichotomy to the one we observed in the case when  In other words, enabling a neutral mutant to live on a more densely connected lattice has a comparable effect on its fixation probability as giving it a certain relative reproductive advantage.Formally, given a population size N and two lattices L N , L 0 N we define the effective fitness, denoted r(L N , L 0 N ), as the unique number r such that In other words, the effective fitness is such a number r(L N , L 0 N ), that a neutral mutant on a lattice L N invading a lattice L 0 N has the same fixation probability as a mutant with relative reproductive advantage r(L N , L 0 N ) in a well-mixed population.
For pairs of low-dimensional lattices with different connectivities d, d 0 , the effective fitness can be computed from the data presented above, see figure 8.We observe that while the effective fitness depends on the connectivities d, d 0 of the two lattices and on their dimensionality, it is mostly independent of the population size N.

Discussion
In this work, we studied the effect of mutations that, rather than altering the reproductive rate of the affected individual, alter how the individual experiences the population structure.To that end, we considered a powerful framework based on the classical Moran birth-death process on graphs, in which the two types of individuals (the novel mutant and the existing residents) perceive the population structure through different graphs.Past evolutionary graph models have discussed inclusion of motility potential as a secondary parameter, besides reproductive rate or fitness [54,56].Our approach generalizes this to the case where the whole migration matrix is genotype dependent, and thus population structure is determined by two graph structures, each for one competing genotype.
As the key quantity, we studied the probability ρ(G A , G B ) that a single neutral mutant which perceives the population structure as a graph G A successfully invades the population of residents which perceive the population structure as a graph G B .For small population sizes, we computed the pairwise fixation probabilities numerically, and we observed that ρ(G A , G B ) tends to be higher when G A is regular and when it contains many edges (that is, the mutant is more motile).We note that the latter aspect contrasts with other models of motility, where an increased dispersal potential of the mutant generally diminishes the fixation probability [55,57,58].
Next, motivated by island models, we considered two regular graphs with the same total number of edges, and we showed that the corresponding fixation probabilities are asymptotically different.In particular, as the population size N increases, the fixation probabilities decay at different rates.Thus, the classic isothermal theorem of [12] does not translate over to the setting in which population structures are type-specific.
Finally, we studied the biologically relevant cases of oneand two-dimensional lattices and we showed that the dispersal radius has similar effect on the fixation probability as the reproductive rate.Recall that in large unstructured populations, a beneficial mutation fixates with constant probability, whereas the fixation probability of a deleterious mutation is exponentially small.Likewise, neutral mutants on lattices with larger dispersal radius have a constant chance of successfully fixating, whereas having lower dispersal radius leads to fixation of the mutant only with exponentially small probability.Thus, in terms of the fixation probability of the mutant, perceiving the population through a more densely connected lattice is effectively equivalent to having an increased reproductive rate.
We view the presented model as a generalization of selection dynamics that is not only needed, but also necessary.For example, for the multicellular structures organized in tissues, the extracellular matrix is produced by a complex mechanism which involves the epithelial cells themselves.To highlight the effect of structural differences between two competing types, we primarily considered the case that mutants might not have different reproductive rate, but have a somewhat altered extracellular matrix.An important example is evolution of carcinoma in epithelial tissues.Malignant tumour cells cannot form the same ligands needed for ECM in a particular tissue.Instead, the organization of tumour cells is royalsocietypublishing.org/journal/rsifJ. R. Soc.Interface 20: 20230355 often reminiscent of random aggregation of cells or an unstructured mesh.This difference in population structure of normal cells and tumour cells can confer significant advantage, or disadvantage, for the tumour cells invading the tissue.The tissue mechanics has been studied in physical models in the past.Recently, selection dynamics has been studied using biomechanical models of tissue formation [68].
In this work, we consider the birth-death dynamics where an individual is selected for reproduction and the offspring replaces one neighbour.Alternatively, one could consider a death-birth version of the model, where an individual is chosen to die first and then one neighbouring cell replaces the void created with an offspring [12,18].Notice that even in the neutral limit of r A = r B , the two update rules would yield different results [69][70][71][72].This has to do with the fact that effective fitness of each genotype depends not only on the reproduction rate, but also on the number of neighbours each genotype sees in its neighbourhood.On average, if a type A perceives more neighbours as a given location it will have higher chance to replace a neighbouring type B. We leave the study of deathbirth process in this context for future work.
Moving on to more complex (though perhaps less realistic) population structures, many natural questions arise.We conclude by commenting on three of them.Recall that for any graph G N on N nodes we have ρ(G N , G N ) = 1/N [61].
First, figure 4 suggests that rðG A N , K N Þ , 1=N for all mutant graphs G A N = K N .While we can prove that ρ(M N , K N ) < 1/N for a graph M N that misses a single edge (see theorem A.1 in appendix A), the general claim is left as an open problem.Similarly, we do not know whether rðK N , G B N Þ .1=N holds for all resident graphs G B N = K N (we do know that it holds for G B N ¼ M N ).Second, following the game theory perspective, Melissourgos et al. [44] asked what is the best mutant response to a given resident graph.That is, given a resident graph G B N on N nodes, which mutant graph G A N on N nodes maximizes the fixation probability rðG A N , G B N Þ?Our results for small graphs show that although the Complete graph K N is frequently the best mutant response, it is not always the case, see figure 9.In particular, when the residents live on a Star graph S 6 , the population is easier to invade through a graph M 6 that misses a single edge, rather than through the Complete graph K 6 -direct computation gives ρ(M 6 , S 6 ) > 0.643 > 0.641 > ρ(K 6 , S 6 ).We note that the difference is minor-both mutant graphs M 6 and K 6 provide a fixation probability well over the neutral threshold value 1/6 ≈ 0.167.
For the complementary question of what is the best resident response G B N to a given mutant graph G A N , the situation is analogous: while the Complete graph is generally hard to invade, it is sometimes not the hardest one.As an example (see figure 9b), when mutants live on a graph G 6 then for the graph G 0 6 we have ρ(G 6 , G 0 6 ) < 0.025 < 0.026 < ρ(G 6 , K 6 ).Third, consider the case when mutants and residents live on different graphs G A N = G B N with comparable edge densities.In this case, we expect that fixation probabilities rðG A N , G B N Þ and rðG B N , G A N Þ typically drop below 1/N.The intuition is that the mutant subpopulation tends to form clusters in G A N but not necessarily in G B N .As a consequence, mutants block each other from spawning onto residents but they do not guard each other from being replaced by residents.However, sometimes the opposite inequality rðG A N , G B N Þ .1=N holds and in general the strength of the effect appears difficult to quantify.In particular, we do not know the answer to the following question: do there exist two regular graphs such that both ρ(G N , G 0 N ) > 1/N and ρ(G 0 N , G N ) > 1/N?We note that if the word 'regular' is dropped, the answer is 'no', as witnessed by the two Star graphs depicted in figure 3.
In conclusion, we have introduced a novel model of stochastic selection dynamics in structured populations, where motility potential for each type is reflected by how that type sees the neighbourhood.This has potential importance for modelling dynamics of aggressive solid tumours in epithelial tissues.The generalized theoretical framework corresponds to a migration matrix or evolutionary graph structure specific to each genotype.Our main observation is that, as a rule of thumb, if a type perceives a higher number of neighbours then it has a selection advantage.We also highlight several counter-examples, especially in small population sizes.Another important observation is that if the two graphs are regular and have the same degree connectivity, the fixation probability is not necessarily 1/N, but its value depends on how the two graphs are structured and overlaid.This might be somewhat surprising, from the perspective of the isothermal theorem.For a more applied case of large lattices, we report the equivalent fitness advantage that a mutant type gains by perceiving the population structure through a more densely connected lattice.Our current model and results can be the basis of a more biophysically inspired model of cancer progression where differences between normal and tumour extra-cellular matrix are included.
Ethics.This work did not require ethical approval from a human subject or animal welfare committee.
Data accessibility.The datasets generated during and/or analysed during the current study are available as part of the electronic supplementary material: http://doi.org/10.6084/m9.figshare.16910170[73].
Declaration of AI use.We have not used AI-assisted technologies in creating this article.

Appendix A. The Complete graph is locally optimal
Here we show that when it comes to ρ(G N ) and r w ðG N Þ, the Complete graph is locally optimal.That is, we show that the graph M N obtained from K N by removing a single edge satisfies ρ(M N ) = ρ(M N , K N ) < 1/N and r w ðM N Þ ¼ rðK N , M N Þ .1=N.
Theorem A.1.Fix N ≥ 2 and let M N be a graph obtained from the Complete graph K N by removing a single edge.Then We claim that wða, bÞ is a solution to this system.Clearly, the formula satisfies w(0, 0) = 0 and w(n − 2, 2) = 1.Since the equation in the last item in the list above can be rewritten as it suffices to check that given an arbitrary configuration (that is, any 0 ≤ a ≤ n − 2 and 0 ≤ b ≤ 2), the expected value of the formula after a single transition does not change.
Ignoring the shared denominator (n − 1) 2 , the expected change in the value produced by the formula in the respective cases is as follows: Grouping the terms with n raised to the same power, the terms in the square brackets on the right-hand side can be rearranged as and finally Since Z = 0 for any admissible integer b (recall 0 ≤ b ≤ 2), we are done.Regarding r w ðM N Þ, a completely analogous proof establishes a formula which implies the desired Since 4k(n − k) ≤ n 2 for any k ¼ 0, . .., n (and with equality only for k = n/2), we can bound On the other hand, the probability p − (k, n − k) that in the next time step we lose a mutant equals Since this expression is independent of k and n, plugging it in the standard formula for the absorption probability of a onedimensional Markov chain we get where e 8 2:718 is Euler's number, we thus obtain Regarding r w ðB N Þ, a completely analogous proof yields  Proof.By a configuration, we mean a subset of vertices occupied by mutants.We denote the possible configurations as a sequence of numbers (corresponding to blocks of consecutive mutants) and symbols ' °' (corresponding to individual residents).That is, for instance the notation k °1 denotes a configuration with k consecutive mutants, then one resident, then one more mutant (and residents before and after).
Proof of the upper bound.This is straightforward.We say that a step of the Moran process is active if it changes the configuration.Given a single mutant, there are two possible active steps leading to immediate extinction (each occurs with probability equal to 1/N • 1/2).On the other hand, there are four possible active steps where mutants reproduce (each occurs with probability equal to 1/N • 1/4).Thus, in total, with probability (2/2)/(2/2 + 4/4) = 1/2 the first active step results in the mutant extinction, hence r 1D N ð4, 2Þ 1=2.
Accounting for trajectories that never reach a configuration with more than m mutants, we can push this upper bound lower: for instance, taking m = 2, denote by x, y, z the extinction probabilities from configurations 1, 2, 1 °1, respectively.Suppose that N ≥ 5 and we are currently at configuration 1.With probability (1/N) • 2/2, one neighbour of the single mutant is selected for reproduction and the offspring replaces the mutant, resulting in mutant extinction.With probability (1/N) • 2/4 the mutant reproduces onto a neighbour and we reach a configuration 2. Similarly, with probability (1/N) • 2/4 the mutant offspring 'skips' the neighbour, resulting in configuration 1 °1.Otherwise, we remain at configuration 1.Thus, conditioning on the first active step (figure 10a) we obtain For k ≥ 2 we have (figure 10b) For analogous expressions of b k , we would need to introduce new variables m k , x k , y k , z k to describe fixation probabilities from configurations listed below the dashed line.Instead of computing the fixation probability exactly, we bound it by considering a different processes M ↓ , where mutants have lower fixation probability than in the original process M. We define M ↓ as follows: the two processes coincide, except that when M reaches any configuration below the dashed line, we remove several mutants so as to obtain a configuration a i or b i (for some i), see figure 10b.For instance, when the current configuration is k °1, the resident from the gap is spawning an offspring, and the offspring moves one place left resulting in a configuration (k − 1) °°1, we additionally remove the single separated mutant, reaching the configuration k − 1, where the mutants have fixation probability a k−1 .Since for any two configurations C ⊆ C 0 we have ρ(C) ≤ ρ(C 0 ), the fixation probability in M ↓ is indeed decreased, as compared with M. (We note that by considering a process M ↑ where we occasionally add mutants so as to only visit configurations a i , b i , one can obtain a stronger upper bound than the one presented above.) In M ↓ we thus obtain It remains to compute the fixation probability in M ↓ , in the limit N → ∞.We do this by a standard argument.Consider a 'potential function' w which assigns a positive real number to each configuration, defined by w(k) = α k , w(k °1) = c • α k , where α, c > 0 are positive real numbers (to be defined later).By solving a system of two equations we find values a 8 0:860, c 8 0:954 (the roots of certain degree-3 polynomials) for which the (expected) potential does not change in one step of the process (that is, the function w is a martingale), except when at configuration 1.Now for any initial configuration x ≠ 1, run the process starting from x until it reaches either a configuration 1 or the configuration N, and let p x be the probability that the former happens.Then wðxÞ ¼ E½wðxÞ j when the process ends ¼ p x Á wð1Þ þ ð1 À p x Þ Á wðNÞ, which rewrites as p x ¼ wðxÞ À wðNÞ wð1Þ À wðNÞ : Since α < 1 we have lim N→∞ w(N) = 0, thus for x ∈ {2, 1 °1} we can write lim Now using the expressions and plugging them into We consider a process M ↓ where, any time a configuration below the dashed line is reached, we remove a mutant to instead reach one of the configurations above the dashed line.We then compute the fixation probability in M ↓ , which is a lower bound.
weights can be any non-negative real numbers, the space of all edge-weighted graphs is immense.To illustrate the difficulties when considering edge-weighted graphs, here we consider a special case, where N = 3, the residents live on a complete (unweighted) graph K 3 , and the mutants live on an edge-weighted graph G x,y,z As in figure 4, each dot corresponds to a graph G and the orange dots correspond to regular graphs.The shape of the point cloud is roughly preserved even as r increases.

Figure 1 .
Figure 1.In epithelial tissues, different cell types align along different lattice-like structures.

Figure 2 .
Figure 2. Moran process with type-dependent dispersal patterns.In each discrete time step, a random individual reproduces and the offspring proliferates to a neighbouring node.Type A offspring (mutant, blue) migrate along the edges of the blue graph G A , whereas type B offspring (residents, red) migrate along the red edges of G B .The key quantity is the fixation probability ρ(G A , G B ) that a single initial mutant successfully invades the population of residents.

Figure 4 .
Figure 4. Small populations N = 6.The fixation probabilities (a) ρ(G) = ρ(G, K N ) (residents on a complete graph) and (b) r w ðGÞ ¼ rðK N , GÞ (mutants on a complete graph) for all 112 graphs G on N = 6 nodes.Each dot corresponds to a graph G; the orange dots correspond to regular graphs.When G = K N , both ρ(G) and r w ðGÞ are equal to 1/N.Other graphs G 6 on six nodes satisfy ρ(G 6 ) < 1/6 and r w ðG 6 Þ .1=6.

Figure 5 .
Figure 5. Dense regular graphs.(a) In a (complete) Bipartite graph B N and a Two-clique graph T N , each node is connected to N/2 other nodes (here N is even).(b) When the mutant lives on B N , the fixation probability satisfies ρ(B N ) ≈ 0.82 • 1/N.By contrast, when the mutant lives on T N , the fixation probability ρ(T N ) tends to zero faster than 1/N.(c) When the residents live on B N or T N , we have r w ðB N Þ % 1:1 Á 1=N and r w ðT N Þ % 1:4 Á 1=N.

Figure 6 .
Figure 6.Overlaying one-dimensional lattices with different connectivities.(a) A circulation graph Cir d N is a one-dimensional lattice with periodic boundary and connectivity (degree) d.We consider d ∈ {2, 4, 6}.(b) When mutants live on a less connected graph (d 1 < d 2 ), their fixation probability decays to 0 at an exponential rate as N → ∞ (here y-axis is log-scale).(c) By contrast, when mutants live on a more densely connected graph (d 1 > d 2 ), their fixation probability tends to a constant.In both panels, the black dashed line shows the neutral baseline 1/N.The values for N ≤ 13 are computed by numerically solving a large system of linear equations.The values for N ≥ 14 are obtained by simulating the process 10 5 times and reporting the proportion of the runs that terminated with the mutant fixating rather than going extinct.

Figure 7 .
Figure 7. Overlaying two-dimensional lattices with different connectivities.(a) We consider two-dimensional lattices with degree 4 (Vonn Neumann neighbourhood), 6 (triangular grid) and 8 (Moore neighbourhood), and with dimensions a × a and a × (a + 1) for a ¼ 3, 4, . . ., 30.(b,c) Similar to the one-dimensional case, the fixation probability decays to 0 exponentially quickly when d 1 < d 2 , whereas it tends to a positive constant when d 1 > d 2 .The black dashed line shows the baseline 1/N.The values are obtained by simulating the process (at least 10 5 repetitions per data point).
First we focus on ρ(M N ) = ρ(M N , K N ).Denote by w(a, b) the fixation probability starting from a configuration with a mutants among the N − 2 fully connected vertices and b mutants among the other two vertices that miss one edge.The values {w(a, b)| 0 ≤ a ≤ N − 2, 0 ≤ b ≤ 2} are the unique solution to the following system of linear equations:-w(0, 0) = 0; -w(N − 2, 2) = 1; and -wða, bÞ ¼P ða 0 ,b 0 Þ p ða,bÞ!ða 0 ,b 0 Þ Á wða 0 , b 0 Þ, where p (a,b)→(a 0 ,b 0) is the probability that, in a single step of the Moran process, we transition to a configuration with a 0 mutants among the N − 2 fully connected vertices, and b mutants among the other two vertices.

Figure 6
Figure6suggests that when d 1 > d 2 then the expression r 1D N ðd 1 , d 2 Þ remains bounded away from 0 as N → ∞, Below we prove this in the special case (d 1 , d 2 ) = (4, 2), that is, when the mutants can place their offspring on the next two vertices and on the preceding two vertices, whereas the residents can place their offspring only on the next one vertex and on the preceding one vertex.

1 (Figure 10 .
Figure 10.Proof of theorem C.1.(a) Upper bound.Boxes represent configurations.Mutants are shown as blue discs, residents as red crosses.Blue (red) transitions correspond to mutants (residents) reproducing.A fraction a/b denotes that there exist a relevant edges, each pointing from a vertex with degree b.With a constant probability, the random evolutionary trajectory leads to mutant extinction without ever reaching a configuration with three or more mutants.(b) Lower bound.We consider a process M ↓ where, any time a configuration below the dashed line is reached, we remove a mutant to instead reach one of the configurations above the dashed line.We then compute the fixation probability in M ↓ , which is a lower bound.

Figure 11 .
Figure 11.Small populations N = 6 with r ≥ 1.The fixation probabilities ρ(G) = ρ(G, K N ) (top row) and r w ðGÞ ¼ rðK N , GÞ (bottom row) for all 112 graphs G on N = 6 nodes, with r ∈ {1, 1.1, 1.5, 2} (columns).As in figure4, each dot corresponds to a graph G and the orange dots correspond to regular graphs.The shape of the point cloud is roughly preserved even as r increases.
Bipartite graph B N is within a constant factor of K N Recall that for N even we denote by B N the (complete) Bipartite graph with parts of sizes N/2 and N/2.We prove that for N large, both ρ(B N ) and r w ðB N Þ are within a constant factor of 1/N.Recall that numerical experiments in the main text suggest that ρ(B N ) ≈ 0.82/N and r w ðB N Þ % 1:1=N, as N → ∞.in the other part.The probability p + (k, n − k) that in the next time step we gain a mutant equals Proof.First we bound ρ(B N ).Consider a time point at which there are n mutants in total: k of them in one part and n − k